!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! Linear solver wrapper
!!
!! \n
!!
!! This module bundles all the linear solvers, thus simplifying using
!! them and calling the respective constructors. See also the
!! comments in \c mod_linsolver_base.
!!
!! \section general General layout
!!
!! The main criteria used in this module are as follows.
!! <ul>
!!  <li> The solution of a linear system always involves three steps:
!!   <ol>
!!    <li> factor the linear system: \c factor; this step can be
!!    further subdivided into two phases (see later for the details)
!!    <ul>
!!     <li> <tt>factor(phase = "analysis")</tt>
!!     <li> <tt>factor(phase = "factorization")</tt>
!!    </ul>
!!    <li> solve the linear system: \c solve; this function can be
!!    called any number of times following a call to \c factor
!!    <li> clean up: \c clean; this function can be called any number
!!    of times, regardless of any previous call to \c factor and \c
!!    solve.
!!   </ol>
!!   These steps are specified as deffered type bound procedures of
!!   the basic type \c c_linpb, so that each solver will specify them.
!!
!!  Concerning the factorization step: the two phases are introduced
!!  because sometimes one has to solve different linear system with the
!!  same sparsity pattern but different matrix coefficients. In this
!!  case, it is convenient to perform the analysis of the symbolic
!!  matrix once, followed by a factorization every time the
!!  coefficients change. The complete interface of \c factor is thus as
!!  follows:
!!  <ul>
!!   <li> if the optional argument \c phase is not present, both
!!   analysis and numeric factorization are executed
!!   <li> if the optional argument \c phase is present, then it must be
!!   either <tt>"analysis"</tt> or <tt>"factorization"</tt>
!!   <li> each call with <tt>phase = "factorization"</tt> must be
!!   preceded by an analysis phase
!!   <li> once an analysis has been performed, there can be an
!!   arbitrary number of calls with <tt>phase = "factorization"</tt>
!!   <li> once an analysis has been performed, any subsequent call
!!   which requires an analysis must be preceded by a call to \c
!!   clean.
!!  </ul>
!!  <li> To define a linear solver, it is necessary to provide an
!!  abstract type extending \c t_linpb; examples are \c c_itpb and \c
!!  c_mumpspb. This type should implement the three subroutines \c
!!  factor, \c solve and \c clean and add all the required fields
!!  which will allow the final user to define his/her linear problem.
!!  <li> The final user will have to define a type which extends the
!!  chosen linear solver, where he/she will provide the required
!!  information. A variable of this type will then be created to
!!  collect the unknown, the linear problem and the solver.
!!  <li> When using MPI based parallel solvers, <em>all</em> the calls
!!  to functions in the present module must be made by all the MPI
!!  processes (but the behavior depends on the chosen solver).
!! </ul>
!! \note Sometimes, the fields provided by \c c_stv are not sufficient
!! for the solver implementation, and one would like to (give the user
!! the possibility to) append information to each \c c_stv object. To
!! avoid introducing an additional type, the solver type derived from
!! \c c_linpb could include an additional field of type \c c_stv so
!! that it could be used as a wrapper of \c c_stv objects in the
!! solver. This idea is used, for instance, in \c mod_newton.
!!
!! \section solvers Supported solvers
!!
!! Currently the following options are implemented (refer to the
!! specific modules for further details):
!! <ul>
!!  <li> <a
!!  href="http://www.cise.ufl.edu/research/sparse/umfpack/">UMFPACK</a>
!!  solvers (direct)
!!  <li> <a href="http://graal.ens-lyon.fr/MUMPS/">MUMPS</a> solvers
!!  (direct), <em>requires MPI</em>
!!  <li> <a href="http://pastix.gforge.inria.fr/">PaStiX</a> solvers
!!  (direct), <em>requires MPI</em>
!!  <li> the iterative solvers of <tt>mod_iterativesolvers</tt>.
!! </ul>
!!
!! \subsection ddc_solvers Parallel execution
!!
!! Some solvers support a parallel solution of a given linear system,
!! other do not. In general, when called by multiple MPI threads, the
!! behaviour of the solver routines will depend on the chosen solver.
!! Currently we have that
!! <ul>
!!  <li> UMFPACK: each thread solves locally the specified system and
!!  there is no communication.
!!  <li> MUMPS: all the threads solve in parallel a unique linear
!!  system; the system matrix cab be centralized or distributed, as
!!  specified in \c linsolver_mumps_parms and documented in the MUMPS
!!  manuals.
!!  <li> PaStiX: all the threads solve in parallel a unique linear
!!  system, the matrix is distributed.
!!  <li> The iterative solvers of <tt>mod_iterativesolvers</tt>: all
!!  the threads solve in parallel a unique linear system.
!! </ul>
!! \warning While in serial execution all the solvers should deliver a
!! similar solution, when using multiple MPI threads a change in the
!! linear solver might imply a change in the problem solved.
!!
!! \section usage Using a solver
!!
!! The general method to define a linear problem with an associated
!! solver is to define a type which extends the chosen solver abstract
!! type and to declare a varable of the extended type. When more than
!! one linear problem must be defined, they can use different types
!! (if the solvers are different) or different variables of the same
!! type (if the linear problems are different but the solver is the
!! same).
!!
!! To be able to select the linear solver during the execution, one
!! will typically
!! <ul>
!!  <li> define a type for each solver that must be available (all
!!  these types will extend, indirectly, \c c_linpb)
!!  <li> define an <tt>allocatable</tt> variable of type \c c_linpb
!!  which will then be allocated to the type of the chosen solver.
!! </ul>
!! The main point here is that, once the variable has been allocated,
!! the calls to \c factor, \c solve and \c clean do not depend on the
!! chosen solver.
!!
!! Here is an example:
!! \code
!! ! Use the relevant modules
!! use mod_state_vars, only: &
!!   c_stv
!! use mod_linsolver, only: &
!!   c_linpb, c_itpb, c_mumpspb
!!
!! ! Define a type for the system unknown
!! type, extends(c_stv) :: t_x
!! contains
!!  ! x_incr, x_tims, x_copy must be implemented
!!  procedure, pass(x) :: incr => x_incr
!!  procedure, pass(x) :: tims => x_tims
!!  procedure, pass(z) :: copy => x_copy
!! end type t_x
!!
!! ! Define the linear problems according to the requirements of the
!! ! available solvers (in practice, provide the abstract procedures)
!! type, extends(c_itpb) :: t_linpb_it
!! contains
!!  ! it_pres, it_pkry, it_scal must be implemented
!!  procedure, nopass :: pres => it_pres
!!  procedure, nopass :: pkry => it_pkry
!!  procedure, nopass :: scal => it_scal
!! end type t_linpb_it
!! type, extends(c_mumpspb) :: t_linpb_mumps
!! contains
!!  ! mumps_xassign must be implemented
!!  procedure, pass(s) :: xassign => mumps_xassign
!! end type t_linpb_mumps
!!
!! ! Define the unknown and the linear problem
!! type(t_x) :: x
!! class(c_linpb), allocatable :: linpb
!!
!! ! Initialize the linear solver
!! select case(linsolver)
!!  case('iterative')
!!   allocate(t_linpb_it::linpb)
!!   select type(linpb); type is(t_linpb_it) ! obvious but necessary
!!     ! set all the fields of a t_linpb_it object
!!   end select
!!  case('mumps')
!!   allocate(t_linpb_mumps::linpb)
!!   select type(linpb); type is(t_linpb_mumps) ! obvious but necessary
!!     ! set all the fields of a t_linpb_mumps object
!!   end select
!!  case default
!!   ! error: unknown linear solver
!! end select
!!
!! ! Solve a linear system
!! call linpb%factor('analysis')
!! call linpb%factor('factorization')
!! call linpb%solve(x)
!! call linpb%clean()
!! \endcode
!!
!! \section details Implementation details
!!
!! The linear solver might use some pointer fields to simplify sharing
!! memory consuming data among linear problems, however the pointed
!! variables are not changed by the solvers. All the working arrays
!! are local to the solvers.
!<----------------------------------------------------------------------
module mod_linsolver

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_linsolver_base, only: &
   mod_linsolver_base_constructor, &
   mod_linsolver_base_destructor,  &
   mod_linsolver_base_initialized, &
   c_linpb

 use mod_iterativesolvers_base, only: &
   mod_iterativesolvers_base_constructor, &
   mod_iterativesolvers_base_destructor,  &
   mod_iterativesolvers_base_initialized, &
   c_itpb

 use mod_gmres, only: &
   mod_gmres_constructor, &
   mod_gmres_destructor,  &
   mod_gmres_initialized, &
   gmres

 use mod_mumpsintf, only: &
   mod_mumpsintf_constructor, &
   mod_mumpsintf_destructor,  &
   mod_mumpsintf_initialized, &
   c_mumpspb

 use mod_umfintf, only: &
   mod_umfintf_constructor, &
   mod_umfintf_destructor,  &
   mod_umfintf_initialized, &
   c_umfpackpb

 use mod_pastixintf, only: &
   mod_pastixintf_constructor, &
   mod_pastixintf_destructor,  &
   mod_pastixintf_initialized, &
   c_pastixpb

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_linsolver_constructor, &
   mod_linsolver_destructor,  &
   mod_linsolver_initialized, &
   c_linpb, c_itpb, c_mumpspb, c_umfpackpb, c_pastixpb, &
   gmres

 private

!-----------------------------------------------------------------------

! Module variables

 logical, protected ::               &
   mod_linsolver_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_linsolver'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_linsolver_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if(mod_messages_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if( mod_linsolver_base_initialized .or. &
mod_iterativesolvers_base_initialized .or. &
                mod_gmres_initialized .or. &
            mod_mumpsintf_initialized .or. &
              mod_umfintf_initialized ) then
     call error(this_sub_name,this_mod_name,                     &
       'The specific linear solve modules are already initialize')
   endif
   if(mod_linsolver_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   call        mod_linsolver_base_constructor()
   call mod_iterativesolvers_base_constructor()
   call                 mod_gmres_constructor()
   call             mod_mumpsintf_constructor()
   call               mod_umfintf_constructor()
   call            mod_pastixintf_constructor()

   mod_linsolver_initialized = .true.
 end subroutine mod_linsolver_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_linsolver_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_linsolver_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   call            mod_pastixintf_destructor()
   call               mod_umfintf_destructor()
   call             mod_mumpsintf_destructor()
   call                 mod_gmres_destructor()
   call mod_iterativesolvers_base_destructor()
   call        mod_linsolver_base_destructor()

   mod_linsolver_initialized = .false.
 end subroutine mod_linsolver_destructor

!-----------------------------------------------------------------------
 
end module mod_linsolver

